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Abstract. We study high-resolution hydrodynamic simulations of Milky Way type galaxies 
obtained within the “Evolution and Assembly of GaLaxies and their Environments” (EA¬ 
GLE) project, and identify those that best satisfy observational constraints on the Milky 
Way total stellar mass, rotation curve, and galaxy shape. Contrary to mock galaxies selected 
on the basis of their total virial mass, the Milky Way analogues so identified consistently 
exhibit very similar dark matter profiles inside the solar circle, therefore enabling more accu¬ 
rate predictions for indirect dark matter searches. We find in particular that high resolution 
simulated haloes satisfying observational constraints exhibit, within the inner few kilopar- 
secs, dark matter profiles shallower than those required to explain the so-called Fermi GeV 
excess via dark matter annihilation. 
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1 Introduction 

Discovering the nature of dark matter (DM) is one of the main challenges for physics today. 
We have evidence of the DM gravitational interactions at different scales, from galactic up 
to cosmological scales (see e.g. [1-4]), but its nature remains a mystery. Weakly Interacting 
Massive Particles (WIMPs) are among the most well-motivated DM particle candidates, 
and they are currently searched for with three detection strategies: (a) direct detection, 
based on the measurement of the recoil energy of nuclei hit by DM particles in underground 
experiments, (b) indirect detection, through the search for secondary particles produced in 
the annihilation or decay of DM, and (c) the search for new particles beyond the standard 
model of particle physics at accelerators, in particular at the Large Hadron Collider. 

The prospects for discovering DM particles with direct and indirect searches strongly 
depend on the distribution of DM in the Milky Way (MW). Recently an extensive compilation 
of measurements of the Milky Way rotation curve - i.e. of the circular velocity of astrophysical 
tracers as a function of the distance from the Galactic centre (GC) - has confirmed the 
existence of large amounts of DM within the solar circle [5], and enabled researchers to 
estimate the DM spatial density profile with parametric [6] and non-parametric methods [7, 
8 ]. 

The DM density prohle, p{R) (where R is the distance from the GC), is a crucial 
ingredient for predicting the intensity flux and anisotropy of gamma rays produced by the 
annihilation of WIMPs in the halo of the MW, since the annihilation rate depends on the 
square of the DM particle number density. While the DM density profile at galactocentric 
distances R > 5 — 6 kpc is relatively well constrained by the analysis of kinematical data of 
the MW rotation curve, the DM profile in the inner few kiloparsecs from the GC is subject 
to large uncertainties [7]. In the absence of observational constraints, the most profitable 
recourse is examination of N-body, and, more recently, hydrodynamic simulations. Pure DM 
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N-body simulations predict a DM density profile that behaves like where 7 « 1 in the 
few inner kiloparsecs, as encoded in the so-called Navarro-Frenk-White (NFW) profile [9]. 
Baryonic processes can either lower or increase the DM density in the inner halo [10-16]. 
However, the size of those effects is still a matter of debate. 

In a recent exciting development for indirect dark matter searches, an unexplained excess 
of gamma rays collected with the Large Area Telescope (LAT) - aboard the Fermi satellite 
- from the centre of our Galaxy has been discovered over and above the standard adopted 
astrophysical background [17-26]. In particular, ref. [25] re-assessed the spectral and mor¬ 
phological properties of the GeV excess, taking into account background model systematics 
associated with the Galactic diffuse emission modelling. The striking similarity of the ob¬ 
served gamma-ray excess with the signal predicted from the annihilation of DM particles in 
the halo of the MW makes the DM interpretation very appealing (see e.g. [23, 24, 27, 28]), 
although other viable astrophysical explanations have been put forward [29-37]. 

Among the alternative GeV excess sources, the possibility that the excess originates from 
a series of leptonic outbursts which occurred ~ 0.1 — 1 yr ago has recently been demonstrated 
to be a viable scenario but a quite unlikely one [36] given the set of parameters needed to 
fully account for the spectral and morphological properties of the GeV excess emission. On 
the other hand, the vigorously debated interpretation in terms of the unresolved emission 
from very dim sources, such as for example pulsars [33] and milli-second pulsars [29-32], has 
been very recently corroborated by two independent works [38, 39]. While it is clear that the 
disc population of unresolved pulsars and milli-second pulsars cannot contribute more than 
10% to the excess emission [40], the contribution from a new population associated with the 
Galactic bulge seems instead to be sufficient to explain the full signal [38, 39]. 

Given the lively debate on possible explanations, and the difficulty of firmly confirming 
or refuting them, it is crucial to fully exploit state of the art simulations to examine what 
is the expected DM profile of MW-like galaxies, in order to refine predictions for the DM 
annihilation gamma-ray flux. We will therefore compare the predictions from the set of 
selected MW-like galaxies with the GeV excess gamma-ray measurements. 

In this work, we study the distribution of DM in MW-like galaxies simulated within 
the EAGLE project [41, 42] - a suite of cosmological, hydrodynamic simulations calibrated 
to reproduce the observed distribution of stellar masses and sizes of low-redshift galaxies 
and designed to address many outstanding issues in galaxy formation such as metallicities of 
galaxies, properties of the intergalactic medium and the effect of feedback on scales ranging 
from dwarf galaxies up to giant ellipticals. We consider at first all galaxies within haloes 
with virial mass in the range 0(10^^ — 10^^) Mq, and we post-process this sample of MW 
analogues by requiring that they satisfy observational constraints on the Galactic rotation 
curve [5], the total stellar mass and the presence of a dominant disc in the stellar component. 

We then evaluate the DM density profile of the final set of what we define to be MW- 
like galaxies and discuss the prospects for DM indirect detection (the implications for direct 
detection will be discussed in an upcoming paper [43]). In particular, we discuss the im¬ 
plications of the resulting DM density profiles for the DM interpretation of the Fermi GeV 
excess. 

After presenting the set of cosmological hydrodynamic simulations used in section 2, we 
will describe our selection procedure and derive the set of MW-like galaxies in our samples 
in section 3. We will dedicate section 4 to the analysis of the DM density profiles of the set 
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Name 

L (Mpc) 

N 

mg (M 0 ) 

mdm (M 0 ) 

e (pc) 

EAGLE HR 

25 

2 X 752^ 

2.26 X lO'^ 

1.21 X 10® 

350 

EAGLEIR 

100 

2 X 1504^ 

1.81 X 10® 

9.70 X 10® 

700 

APOSTLEIR 

- 

- 

1.3 X 10® 

5.9 X 10® 

308 

APOSTLE HR (I) 

- 

- 

1.0 X 10^ 

5.0 X 10"^ 

134 

APOSTLE HR (H) 

- 

- 

5.0 X 10® 

2.5 X 10"^ 

134 


Table 1. Parameters of the simulations discussed in this paper. L is the comoving sidelength of the 
simulation cube, N the number of simulation particles prior to splitting, Wg the initial gas particle 
mass, rudm the DM particle mass, and e the Plummer-equivalent physical softening length. The 
resolution limit is usually taken to be 2.8xe, i.e. 1.96, 0.98 and 0.87 kpc for EAGLE IR, EAGLE HR 
and APOSTLE IR, respectively. 


of MW-like galaxies and finally section 5 to the discussion of the implications for the GeV 
excess emission. We conclude in section 6. Appendices A and B contain additional material 
supporting our findings. 

2 Simulations 

In this section we briefly describe the set of simulations used, which form part of the EAGLE 
project [41, 42]. The EAGLE simulations were performed using a modified version of the P- 
GADGEtS Tree-Smoothed Particle Hydrodynamics (SPH) code [44] that has been modified 
to use the state-of-the-art SPH flavour of [45] (whose impact on the galaxy population is 
discussed by [46]). The cosmological parameters were chosen to be those derived from the 
analysis of the Planck 2013 measurements [47] and they have the following values: 12^ = 
0.307, r^A = 0.693, rife = 0.0482, h = 0.678, = 0.83, and Ug = 0.961. The simulations 

were run at two different resolutions, which we refer to as intermediate (EAGLE IR) and 
high resolution (EAGLE HR). The former was run in a series of cosmological volumes up to 
a maximum of 100 Mpc on one side, and the latter up to 25 Mpc. The simulations start 
at z=127 with an equal number of gas and DM particles whose masses are given in table 1. 
Gas particles have a finite probability to be turned into star particles that increases with the 
gas pressure, such that the local star formation law [48] is reproduced. Each star particle 
represents a simple stellar population, and inherits the mass and element abundances of the 
parent gas particle. The properties of the two EAGLE runs used in this study (EAGLE IR 
and EAGLE HR) are reproduced from [41] in table 1. In addition to the EAGLE cosmological 
volumes, we also make use of simulations of the APOSTLE project [49, 50]. These simulations 
use the same code as the EAGLE project applied to zoomed regions containing a close pair 
of ~ IO^^Mq dm haloes that could host our MW galaxy and M31, i.e. be an analogue 
of the Local Group. We use twelve APOSTLE volumes simulated at similar resolution to 
EAGLE HR, which as a group we denote APOSTLE IR (for “intermediate” resolution). In 
addition, APOSTLE consists of other two re-simulations at ten times higher mass resolution, 
denoted generically as APOSTLE HR. Their simulation details are also included in table 1. 
One minor difference from the EAGLE cosmological volumes is that they use the WMAP-7 
cosmological parameters: Dm = 0.272, Da = 0.728, D^ = 0.0455, h = 0.704, = 0.81, and 

Ug = 0.967. 

To select MW-like galaxy candidates, we first extract all galaxies located at the minimum 
of their halo potential wells as returned by the Subfind algorithm [51] (i.e. we exclude 
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satellite galaxies) in haloes of virial mass, M 200 , in the range 5 x 10^^ < M 200 /M.Q < 1 x 10^^, 
where M 200 is defined as the mass enclosed within the sphere that contains a mean density 200 
times the critical density. We extend our range to the most massive haloes because the model 
underpredicts slightly the stellar masses contained within haloes of M 200 rsj IO^^Mq [41], 
and therefore MW mass galaxies are found in haloes that are slightly more massive than 
is inferred from abundance matching results, e.g. [52], We show in appendix A that this 
mismatch between halo mass and stellar mass has little effect on our results.^ 

3 Selection of Milky Way-like galaxies 

We consider the simulation runs EAGLE IR, EAGLE HR and APOSTLE IR described in 
section 2. We start from the corresponding subsets of galaxies at the centre of haloes with 
5 X 10^^ < M 200 /MQ < 1 X 10^^. The initial sets are composed of 2411, 61 and 24 objects 
for the EAGLE IR, EAGLE HR and APOSTLE IR run, respectively. 

In this section we aim to select the galaxies that most closely resemble the MW. Our 
definition of MW-like galaxies is based on a minimal set of criteria that the simulated haloes 
should satisfy. In particular, for a simulated halo to host a good MW-like galaxy, we require 
that: 

(i) The simulated rotation curve fits well the observed MW kinematical data in ref. [5]. 
We explain the method followed to derive the rotation curves from the simulation, the 
data used in the analysis and the goodness of fit definition in section 3.1. 

(ii) The total stellar mass of the simulated galaxies is within the 3a MW range derived from 
observations, 4.5 x 10^*^ < M*/Mq < 8.3 x 10^° [53]: 335, 12, and 2 galaxies satisfy this 
constraint in the EAGLE IR, EAGLE HR and APOSTLE IR run respectively.^ 

(hi) The galaxies contain a substantial stellar disc component. See section 3.2. 

The three criteria listed above define our MW-like galaxies. The final sets of good objects 
are presented in section 3 for the EAGLE IR, EAGLE HR and APOSTLE IR simulations. 

We are aware of the fact that besides structural parameters, such as the rotation curve 
used in the present work, more criteria should be imposed to identify truly MW-like galaxies, 
such as brightness profile constraints, star formation history, metallicity gradient, and disc 
scale length/height. However, we stress that our main purpose is to derive implications 
for indirect DM searches rather than to test the ability of the EAGLE simulation code to 
reproduce the MW. Indirect detection prospects depend on the shape of the DM profile, 
which turns out to be almost universal for the simulated objects in EAGLE IR, EAGLE 
HR and APOSTLE IR as shown in ref. [54] and more thoroughly in section 4. The method 
presented here is general and can be applied to future simulations. 

3.1 Observed Milky Way rotation curve and goodness of fit 

Recently, ref. [5] collated a large amount of observational measurements of the MW rotation 
curve and compared these with the expectations from a large set of baryonic models, finding 

^Appendices refer to the EAGLE set of simulations. However, the conclusions reached there can be safely 
extended to the APOSTLE IR run. 

^We note that for APOSTLE IR one of the two galaxies falls on the lower boundary of the 3cr observed 
range. In order to have more than one galaxy in the final selection, we keep this one as well. 
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robust evidence for DM even within the inner 5-6 kpc of the Galaxy. We make use of the 
recent compilation of kinematic tracers presented in [5] . We adopt a local circular velocity of 
vq = 230 km/s, a local galactocentric distance of i?o = 8 kpc, and the component of the solar 
peculiar velocity in the direction of the Galactic rotation, Vq = 12.24 km/s [55], as fiducial 
values for the analysis presented below. However, owing to the thorough discussion by [5] 
about the uncertainty in the rotation curve data due to vq and Rq, we dedicate appendix B 
to the scrutiny of possible variations in the results due to different choices of vq and Rq and 
we show that our main conclusions remain unchanged. 

The observational data are provided as constraints on the angular circular velocity, 
0Jc{R) = Vc{R)/R, and the galactocentric distance R. Here Vc{R) is the circular velocity at 
distance R. As explained by [5], using the angular circular velocity for fitting purposes is 
more convenient than working with Vc{R), since the errors of uJc{R) and R are not correlated, 
while this is not the case for Vc{R) and R. 

The circular velocity, Vc{R), is defined as the velocity of a test particle on a circular 
orbit at radius R from the GC. By equating the centripetal and gravitational forces on the 
test particle, it is simple to show that, for a spherically symmetric matter distribution: 


(3.1) 

where M(< R) is the total mass enclosed within radius R, and G is the universal gravitational 
constant. In figure 1, we display the kinematical data in the plane (uic, R) together with the 
rotation curves predicted by all haloes in the EAGLE IR (top left panel), EAGLE HR (top 
right panel), and APOSTLE IR (bottom panel) simulation runs, satisfying our halo mass 
constraint. Erom figure 1, it is evident that there is a wide range of variation in the rotation 
curves from the simulation when considering all objects whose halo mass, M 20 O) lies in the 
selected MW mass range, 5 x 10^^ < M 2 OO/M 0 < 1 x 10^^. However, by forcing the haloes 
to have the correct total stellar mass (in the 3 ct observed range for the MW), the number 
of good objects is reduced significantly: from 2411 to 335 for EAGLE IR, from 61 to 12 for 
EAGLE HR, and from 24 to 2 for APOSTLE IR. We also display the haloes that give the 
best fit to the kinematical data (see below for further details). Therefore, classifying a halo 
with the correct halo mass as MW-like is too simplistic a criterion, and thus will often fail to 
reproduce the MW kinematical data. For this reason, we extend the definition of MW-like 
galaxies to account for the agreement with the observed Galactic rotation curve.^ 


To derive the goodness of fit of the simulated rotation curve to the observed data, it 
is convenient to work with the reduced quantities x = R/Rq and y = uJcjuiQ — 1 [5], where 
Wo = vq/Rq. The two-variable can be written as: 


2 ^ iVi - y{xi)f 

^ + {dyi/dxiY al. ’ 


(3.2) 


^More realistically the requirement of spherical symmetry is not guaranteed to be a good approximation, 
and thus Equation 3.1 may not be appropriate. As a check, for the APOSTLE IR galaxies we calculated a 
series of rotation curves by summing up the gravitational forces due to all of the simulation particles in the 
box. We found that the difference in rotation curve amplitude between the full calculation and that obtained 
from the spherical symmetry is typically less than 10% at all radii. Moreover, we checked that the values 
of these new rotation curves (for the ht to the measured MW rotation curve) are not signihcantly different 
from the ones obtained by assuming spherical symmetry. Most importantly, the ordering of haloes based on 
values stays the same. We therefore conclude that our assumption of spherical symmetry is appropriate 
for this study. 
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Figure 1 . Kinematical data from [5] used in this work {black points and error bars) and rotation 
curves of simulated haloes for the EAGLE IR {top left panel), EAGLE HR {top right panel), and 
APOSTLE IR {bottom centre panel) runs. The green, magenta, and orange curves correspond to 
galaxies which fit our criteria as defined in section 3 and give the best fit (lowest in eq. (3.2)) in the 
EAGLE IR, EAGLE HR, and APOSTLE IR runs, respectively. The light and dark coloured bands 
correspond to all simulated objects (5 x 10^^ < M200/MQ < 1 x 10^“^) and those which, in addition, 
satisfy the observed MW stellar mass range (4.54 x 10^° < M*/Mq < 8.32 x I0^°), respectively. The 
vertical dashed blue line marks the minimal radius considered in the fitting procedure (see text for 
further details). 


where i runs over the observational data points considered, and y{xi) is the simulated rotation 
curve evaluated at Xi = Ri/Rq. Both experimental errors in x {(Jx^) and y {o^y^) are considered 
and obtained through standard error propagation from the errors in Wc and R. For the fit, 
we consider measurements at R > 2.5 kpc. Indeed, the data derived by [5] assume circular 
orbits for the tracers and this approximation can break at small radii due to the effect of the 
Galactic bulge gravitational potential. Note also that 2.5 kpc is larger than the gravitational 
softening length used in our simulations. 

In figure 2, we show the reduced versus the total stellar mass, with N=2687 the 
total number of observational data points for the three sets of simulations. From the top 
left panel, which refers to EAGLE IR run, it is evident that the global minimum of the 
X^/(Af — 1) distribution naturally occurs within the 3 cj measured range of the MW total 
stellar mass [53]. In other words, a good match of the simulation with the measnred MW 
rotation cnrve is given by galaxies that have the correct MW stellar mass. In general, the 
contribution to rotation curves from stars (which largely dominate over the gas in the total 
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Figure 2. Reduced X^/i^ ~ 1)) versus total stellar mass for EAGLE IR {top left panel), EAGLE 
HR {top right panel) and APOSTLE IR {bottom panel). The coloured dots indicate the whole set 
of 2411, 61 and 24 simulated haloes in EAGLE IR, EAGLE HR and APOSTLE IR; of those, 335, 
12 and 2 galaxies lie in the 3tT measured range {black dashed vertical lines) of the MW total stellar 
mass [53]. The colour-bar shows the distribution of halo mass, M 2 oo- 


baryonic component) is larger than the DM contribution up to i? < 5 kpc. We note that by 
performing the analysis with a distance cut at 5 kpc our results remain unchanged. 

The halo masses of the simulated galaxies with correct stellar mass and in the minimum 
reduced region are larger than expected from observations of the MW, as for example 
-^^ 200 , MW = X 10^^ M© [56]. This result might indicate that the feedback in EAGLE 

simulated haloes in this mass range is slightly too efficient, and thus the stellar mass per 
unit halo is suppressed with respect to that inferred from estimates of the MW stellar and 
halo masses. This shortcoming is reflected in the stellar luminosity function, in that EAGLE 
underpredicts the abundance of galaxies with stellar masses ~ 2 —8x 10^*^M© relative to what 
is observed in galaxy surveys [41]. We also note that the total MW halo mass is affected by 
large uncertainties, with estimates based on kinematics of satellites, abundance matching, and 
the local Hubble flow, yielding somewhat discrepant results [57-60]. However, in appendix A, 
we show that the quantities relevant for DM indirect detection, and, in particular, for the 
implications on the GeV excess interpretation, are not affected by the large halo mass. We 
are thus confident that this mismatch between halo mass and stellar mass has little effect on 
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Figure 3. Distribution of stellar circularities, /(e), for the haloes selected in the EAGLE IR {top 
left panel), EAGLE HR {top right panel), and APOSTLE IR {bottom centre panel) resolution runs 
according to our criteria, including the requirement of having a significant disc stellar component. 
Thick lines refer to the haloes giving the best-fit to the Galactic rotation curve for the three runs. 


our final results. 


3.2 Morphology of simulated galaxies: disc and spheroid 


The MW is a spiral galaxy, with a well defined, unperturbed thin disc, a small bulge and 
a bar component. We therefore seek to select simulated galaxies that are themselves disc 
galaxies rather than ellipticals or undergoing mergers. 

Following the approach of ref. [61], we characterise the dynamics of each simulated 
galaxy by looking for evidence of coherent rotation. Each star particle in the simulated 
galaxies possesses an angular momentum vector relative to its host’s standard of rest. Any 
bulk stellar component that has the same angular momentum vector as that of the hosting 
simulated galaxy will be considered to belong to the disc. We can therefore use the dis¬ 
tribution of angular momentum vectors of individual particles relative to the net angular 
momentum of the galaxy to discriminate between discs (coherent rotation) and spheroids 
(no coherent rotation; comprises bulges and stellar haloes). 

For each selected galaxy, we derive the distribution of the stellar orbital circularity 
parameter, e: 


e(r) = 


Jz 

jc{r) 


Jz 

ruc(r)’ 


(3.3) 


where jz is the component of the specific angular momentum parallel to the total angular 
momentum of the galaxy, and jc(^) is the total specific angular momentum of a circular orbit 
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Name 

M*[xI0^'J Mq] 

M2 oo[xIO^^ Mq] 

D/T 

xV(^-i) 

EAGLE IR 

6.69 

4.38 

0.70 

16.56 


6.95 

2.39 

0.61 

17.32 


7.16 

9.12 

0.69 

17.59 


6.10 

2.61 

0.51 

17.67 


7.91 

6.28 

0.52 

18.23 


5.53 

5.36 

0.50 

18.30 


6.78 

2.33 

0.60 

18.47 


6.22 

2.23 

0.69 

18.83 


5.50 

3.09 

0.65 

18.84 


6.83 

3.08 

0.43 

19.45 

EAGLE HR 

5.31 

3.50 

0.45 

74.55 


5.48 

2.76 

0.46 

220.96 

APOSTLEIR 

4.48 

2.15 

0.50 

51.04 


4.88 

1.68 

0.70 

221.27 


Table 2. Relevant parameters of the finally selected MW-like galaxies that satisfy our selection 
criteria in the EAGLE IR, EAGLE HR and APOSTLE IR runs. We quote: total stellar mass, halo 
mass, disc-to-total mass ratio (D/T), and reduced for the fit to the rotation curve data. 


at distance r; here specific angular momentum is defined as angular momentum divided by 
the star particle mass. The distribution of stellar circularities, /(e), is a good indicator of 
the relative importance of the disc with respect to the spheroidal stellar component. A disc 
in rotational support, that is a configuration in which gravitational collapse is offset by the 
centripetal acceleration, corresponds to a distribution peaked at about e = 1, while stars 
in a system supported by dispersion (i.e. a spheroidal system) show an almost symmetric 
distribution around e = 0 [61-63]. 

As a further constraint - in addition to the goodness of fit to the observed MW rotation 
curve - we want to select objects whose stellar kinematics shows a disc component and, 
thus, are not completely dominated by the spheroidal contribution. When building /(e), 
we calculate the net specific angular momentum of the disc using only those star particles 
with galactocentric radii in the range 3-20 kpc. In this way, our determination of the disc 
angular momentum direction should neither be affected by the isotropic motions of bulge 
stars at low radii nor by high angular momentum halo stars at large radii. We inspect the 
distribution of stellar circularities, /(e), for the haloes giving the best fit to the rotation 
curve data with the aim of building a final snbset of objects that pass our criteria and can 
be considered to be MW-like in our perspective. We retain a galaxy if the stellar fraction in 
the range e > 0.45 is larger than 50%. This criterion is meant to identify galaxies that have 
a dominant disc and, on the other hand, to remove galaxies that show an almost symmetric 
distribution around e = 0 - and can thus be classified to have a spheroidal only component. 
We tested different cuts on /(e) and adopted the one that is most conservative in retaining 
objects with a significant disc component. We note that the above criterion relies on the 
stellar kinematics only, and that it might differ from the bulge-to-disc decomposition given 
by photometric measurements. Indeed, photometric observations depend on assumptions on 
the shape of the brightness profile of the spheroidal and disc components [62, 64]. However, 
a full comparison of the simulation outcome and the photometry of the MW is beyond the 
scope of the present work and it would be interesting to study with follow-up analyses. 
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In the case of the EAGLE IR run the number of galaxies passing the stellar mass and 
disc dominance criteria is 145. Therefore, we further reduce the set of objects of the EAGLE 
IR resolution run by selecting those that have a reduced < 20 for the fit to the MW 
rotation curve data. Eor the EAGLE IR run, we thus finally select 10 MW-like galaxies. Eor 
the EAGLE HR resolution run, instead, we do not impose any further constraint on the 
since only 2 objects survive the /(e) criterion (and lie in the correct stellar mass range), due 
to the 64 times smaller simulation volume. Finally, for the APOSTLE IR run, the 2 galaxies 
satisfying the dynamical and stellar mass constraints also match the disc requirement. 

In figure 3, we show the distribution of stellar fraction for the galaxies that give the best 
to the rotation curve data, have the acceptable MW stellar mass, and show a substantial 
stellar fraction in the disc in the EAGLE IR, EAGLE HR, and APOSTLE IR runs. The 
reduced x^ values are quoted in table 2. We note that in the cases of the EAGLE HR 
and APOSTLE IR runs, the statistics for the initial sample are quite poor. In these cases 
the galaxy with the lowest x^ does not have any privileged statistical meaning, since the 
true minimum of the distribution is not guaranteed to have been found because of the small 
statistics. However, from the EAGLE IR results we are confident that the global minimum 
should lie within the correct MW stellar mass range. The fact that the values of the 
EAGLE HR and APOSTLE IR runs are worse than the EAGLE IR one is thus not surprising. 
In what follows, we will use the best-fit galaxy as a reference for illustrating the implications 
for DM indirect detection. 

For the sake of completeness, we also quote here the reduced x^ of the 4 galaxies 
belonging to the APOSTLE HR which are analysed in detail in ref. [65]. The reduced 
X^ values are (in order): 1203, 1982, 2056, 3554. Those are clearly poorer than the runs 
analysed in this work, but, again, the sample is very limited in number. Moreover, we have 
checked that those galaxies have a lower total stellar mass than what is observed within 3a 
(M*/(10^°Mq) = 2.22, 1.56, 1.05, and 1.06, from smallest to largest x^ respectively). Since 
the stellar mass plays a role in the determination of rotation curves in the inner 5 kpc, such 
a discrepancy in the stellar mass can be a further reason for the poor x^- 

To corroborate the implementation of criterion (iii), in figures 4, 5 and 6 we display 
the images of the edge-on and face-on orientations of the selected MW-like galaxies in the 
EAGLE IR, EAGLE HR and APOSTLE IR run. The images are produced with the radiative 
transfer code skirt [66]. This code generates images of the galaxy in the u, g, and r SDSS 
filters within a 30 kpc aperture of the galaxy centre; the dust distribution is inferred from 
the gas metal distribution in the interstellar medium [67, 68]. From those images one can 
see that the majority of the selected objects shows clearly a disc component. 

To further test the prominence of the disc, we compute the disc-to-total mass ratio for 
the set of selected MW-like galaxies. The disc-to-total mass ratio writes as [61], 


D _ Mg 
T Mrf + M, ’ 


(3.4) 


where Md = Yle >o stellar mass in the disc, and Mg = X]e>-o 6(/(^*) ”^ 9 ) 

is the stellar mass in the spheroid. The disc-to-total mass ratios are listed in table 2 for the 
sets of the selected MW analogues in the three simulations and are all in the range 0.4 - 0.7. 
For the EAGLE HR simulation run, the ratios are systematically lower than what is expected 
for real spiral galaxies [69], while for some galaxies in the EAGLE IR and APOSTLE IR runs 
the D/T ratio is closer to observations of the MW, (D/T)mw ~ 0-86 [53]. 


- 10 - 



In summary, we have identified a minimal set of 10, 2 and 2 galaxies for the EAGLE IR, 
EAGLE HR and APOSTLE IR simulation runs that simultaneously (i) give a good fit to the 
rotation curve of the MW, (ii) have the right total stellar mass (within 3a) of the measured 
value, and (hi) show a significant disc stellar component. In table 2, we summarise the 
main properties of the chosen galaxies for the three simulation runs. These three subsets of 
objects represent MW analogues according to our definition. We have thus demonstrated that 
minimal selection criteria, such as the requirement to give a good fit to the observed rotation 
curve, significantly reduce the number of prototypical MW-like galaxies. In sections 4 and 5 
we will use the two sets of galaxies to investigate the distribution of DM in the halo and 
discuss implications for indirect DM searches. We finally stress that discussing the matching 
of EAGLE simulated galaxies with other MW observables, such as the brightness emission 
profile, is beyond the scope of the present paper. Eurther studies along this direction would 
be extremely helpful for the field. 

4 The Galactic dark matter density profile 

In this section we analyse the DM density profile of the final set of MW-like galaxies and dis¬ 
cuss possible implications for DM indirect detection. As explained in section 1, the predicted 
DM annihilation flux is very sensitive to the DM density profile and it is thus important to 
study the predictions of hydrodynamic simulations and any deviations from pure DM N-body 
simulation density profiles. As for the EAGLE simulations, it has been shown that, on aver¬ 
age, haloes in the MW mass range (~ IO^^Mq) have a DM density profile that significantly 
deviates from the NEW profile and has a core in the inner few kiloparsecs [54]. Although 
the flattening might appear below the effective resolution (depending on the resolution run 
considered), the presence of a core on small radial distances from the GC seems to be a 
universal feature of EAGLE simulated objects. This has been recently demonstrated for the 
set of galaxies of the APOSTLE project, which were run with the same code and subgrid 
physics model as EAGLE; the density profile properties of the highest APOSTLE simulation 
run have been analysed in detail in [65]. We turn now to the analysis of the density prohles 
of our final set of MW-like galaxies. 

In the left panels of figure 7, we show the DM density profile of the final set of haloes for 
the EAGLE IR, EAGLE HR and APOSTLE IR runs. The DM density is derived directly from 
the mass enclosed in a given shell between R and R -|- 6R, adopting equally spaced radial bins 
in logarithmic space. By construction, we assume spherical symmetry for the distribution 
of DM, which has been shown to be a good assumption for the APOSTLE simulations [65]. 
Eurthermore, it has been shown that baryons tend to make the distribution of DM more 
spherical than in simulations including solely DM [70-72]. We verified that for our set of 
MW-like galaxies the assumption of a spherically symmetric profile is a good approximation. 
The uncertainty in the density is given by the Poissonian error in the number of particles 
in each mass shell (the error in the distance refers instead to the adopted binning). As 
presented in section 2, the effective resolution limit of the simulation runs, i.e. the Plummer- 
equivalent softening length, is 2.8 x e, where e = 0.7 kpc, 0.35 kpc and 0.31 kpc for the 
EAGLE IR, EAGLE HR and APOSTLE IR runs, respectively. However, the radius at which 
profiles can be considered as converged is larger than this value and can be estimated in 
collisionless simulations using the criterion of [73] that identifies the radius at which the 
integral in mass is independent on the resolution. The so-called “Power radius” is Rpos = 
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Figure 4 . Pairs of mock three-colour composite grid observations for the 10 simulated galaxies in 
the EAGLE IR run with edge-on (left) and face-on (right) projections, within a 30 kpc aperture of 
the galaxy centre obtained using the skirt [66] radiative transfer code as described by [67, 68]. 


3.6 kpc, 1.8 kpc and 1.8 kpc^ for the EAGLE IR, EAGLE HR and APOSTLE IR runs. The 
very concept of convergence - that ever increasing resolution will cause the measured quantity 
to asymptotically tend towards some finite value - is much less clear in simulations that 

"^The Power radius is a halo-dependent quantity and thus it is defined for each halo. We checked however 
that the variation of Rpos among haloes of the same simulation run varies only within a few percent. Therefore, 
using an average value for Apos is a good approximation. 
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Figure 5. Same as figure 4 for the 2 selected MW-like galaxies in the EAGLE HR run. 



Figure 6. Same as figure 4 for the 2 selected MW-like galaxies in the APOSTLE IR run. 
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Figure 7. DM density profiles {left panels) and the radial change of the local logarithmic slopes 
{right panels) of the selected MW-like galaxies in the EAGLE IR {top), EAGLE HR {middle) and 
APOSTLE IR {bottom) runs. The thick grey line represents the prediction for an NEW profile with 
Tg = 20 kpc and local DM density p© = 0.4 GeV /cm? (as commonly assumed in DM indirect detection 
studies). In all panels the effective resolution of the simulation is shown by the dashed black line, while 
the black arrows on the left panels indicate the convergence radii of 3.6 kpc (EAGLE IR) and 1.8 kpc 
(EAGLE HR and APOSTLE IR) as discussed in the text. 


feature baryon physics, and so the innermost radius at which the profiles may be considered 
converged is ill-defined. A discussion of these issues can be found in refs. [41, 54, 65]. 

In figure 7 we show both the resolution limit and the Power radius for the three resolution 
runs. Between those two radial scales the results of the simulation have to be treated with 
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extreme caution, as there might still be some residual numerical effects due to discreteness 
and softening of the gravitational force. This is particularly true for the EAGLE IR run, 
while it is less dramatic for EAGLE HR and APOSTLE IR. 

In general, the DM density is shallower than what is expected from an NEW profile in 
the inner 1.5 ~ 2 kpc. Between about 1.5 - 6/8 kpc, however, the effect of baryons results 
in a steeping of the DM profile. Analogous features have been found in the APOSTLE HR 
simulations [65]. The origin of the flattening in [65] might be related to the star formation 
history of the haloes: while they developed a cusp at redshifts larger than 1, at ^ < 1 episodes 
of sudden enhancement of the star formation rate occurred. Such violent starbursts could 
have destroyed the cusp [10] at sufficiently late times that the galaxy cannot subsequently 
accrete enough small DM haloes to rebuild the cusp. We stress here that the mild flattening 
occurring in EAGLE IR between the resolution limit and the Power radius is most likely due 
to softening effects and limitations in the gas physics model [65], while for EAGLE HR and 
APOSTLE IR it might be a combination of the softening and of physical effects. 

To better appreciate the deviation of the simulated haloes’ profiles from pure DM sim¬ 
ulations, i.e. from an NEW density profile, in the right panels of figure 7 we show the radial 
variation of the local logarithmic slope for our selected MW-like galaxies together with the 
expectation for the NEW density profile. In the case of the EAGLE IR run, it is not possi¬ 
ble to establish the presence of the flattening since the error bars on the logarithmic slope 
make it compatible with the NEW profile down to the resolution limit. This is compatible 
with what is found by [54], who studied properties of stacked haloes in the EAGLE IR and 
EAGLE HR simulations. However, for the EAGLE HR and APOSTLE IR selected MW-like 
galaxies, there is a deviation of the slope from -1 and a tendency to 0 slightly above the 
resolution limit. From the same panels (EAGLE HR and APOSTLE IR runs) we can also 
appreciate the effect of the baryonic contraction: in the range 1.5 - 6 kpc the logarithmic 
slope is steeper than -1. 

The two features of the DM profile (flattening below 1.5 kpc and adiabatic contractions 
below 10 kpc) appear thus to be universal properties of the simulated MW-like galaxies within 
the EAGLE galaxy formation model, i.e. a generic outcome of the subgrid model assumed. 
What is striking, though not surprising, is that for the selected MW-like galaxies the overall 
normalisation of the different DM haloes does not show large scatter. This is again a direct 
consequence of demanding that the selected galaxies fit the kinematical data of the MW. As 
a consequence, the variation in the local DM density is also small. For the final sets MW-like 
galaxies, p© (with R© = 8.0 kpc) spans the range 0.44 - 0.59 GeV/cm^ (EAGLE IR), 0.41 - 
0.56 GeV/cm3 (EAGLE HR), and 0.41 GeV/cm^ (APOSTLE IR). This is compatible with 
recent measurement of the local DM density, whose standard value at 8.5 kpc is assumed to 
be 0.4 GeV/cm^ (which translates to 0.44 GeV/cm^ at 8.0 kpc) [74, 75]. We acknowledge 
however that the local DM density is affected by several uncertainties and, in particular, the 
DM density has been found to be larger in the stellar disc compared to what is derived from 
a spherical shell [76]. 

5 The implications for dark matter indirect detection 

As explained above, the DM density profile enters squared in the prediction of the gamma-ray 
intensity flux and anisotropy signal that comes from annihilation of DM particles in the halo 
of the MW, see for example [77]. Searches for DM at the GC are thus especially sensitive to 
the uncertainty affecting the determination of the inner distribution of the DM density. In 
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Prohle 

((Tu)[xl 0 ^‘’cm^/s] 

[GeV] 


p-value 

gNFW ( 7 = 1 . 26 ) 

1.71 ±0.11 

47.32 ± 1.07 

223.9 

0.73 

EAGLE HR 

1.96 ±0.14 

46.37 ± 1.37 

246.3 

0.34 

APOSTLEIR 

1.76 ±0.16 

45.36 ±2.96 

283.9 

0.02 


Table 3. Best-fit parameters {{crv) and mass) together with goodness of fit indicators - with 238 
degrees of freedom - and p-value) of the ten sub-regions fit for the gNFW profile, and the best-fit 
galaxies of the EAGLE HR and APOSTLE IR runs. 100% annihilation into b-quarks is assumed. 



[GeV] 


Figure 8 . Constraints on the {{<Jv), m^) plane for annihilation into b-quarks. The three contours 
refer to different DM density prohles: the EAGLE HR run best-fit MW-like galaxy’s DM density 
prohle (magenta), the APOSTLE IR one (orange), and the gNFW profile with 7 = 1.26 (grey) 
(cf. table 3). 


this section, we focus on the implications that the density profiles discussed above entail for 
the DM interpretation of the Fermi GeV excess. 

We remind the reader that the analysis in [25] tested several templates for the GeV ex¬ 
cess spatial distribution, i.e. morphology, assuming an underlying generalised NFW (gNFW) 
profile: 

(r/rs)~'^ (l + r/rs)‘^~'^ ' 

and varying the inner slope 7 in the range 0.7 - 1.5 (thus accounting for both shallower and 
steeper profiles at small radii than the traditional NFW). The best-fit value of 7 was found 
to be 1.26 ± 0.15, corresponding to an inner slope steeper than the standard NFW prohle. 

In order to verify the plausibility of this interpretation, it is crucial to test the prohles 
predicted by hydrodynamic simulations against the GeV excess data. In particular, while 
we know already that the spectral shape of the signal is compatible with a DM annihilation 
signal into light and b-quarks pairs (even when properly taking into account the background 
model systematic uncertainties) [27], the spatial distribution of the signal represents the key- 
point-test for any model under scrutiny. We use the GeV excess data as derived by [25] in 
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the ten segments analysed in the region \l\ <2° and 2° < |6| < 20°. This region corresponds 
to a radial scale from about 0.3 kpc up to 3 kpc. Clearly, the lower limit is not probed by 
any of our three simulations whose effective resolution, 2.8 x e, is at best 0.87 kpc. Therefore, 
to make predictions down to the scale at which the GeV excess is measured we need to 
extrapolate the simulated profiles below the resolution limit. 

Following a traditional approach to the extrapolation issue, we adopt a power-law whose 
steepness is the maximal compatible with the total mass inside the extrapolation radius, 
namely the maximal asymptotic slope [78, 79]. At each radius r, mean and local densities 
define a robust upper limit on the asymptotic inner slope: 7max(?') = 3(1 — p{r)/p{r)). 
However, as already explained in section 4, possible numerical effects still might occur between 
the resolution limit and the Power radius. The Power radius was developed for the specific 
purpose of ascertaining the minimum radius at which the enclosed mass is converged [73], 
and is a useful guide for this test. For this reason, we believe that only extrapolating from the 
Power radius guarantees the results to be stable against residual numerical effects. Moreover, 
we checked that extrapolations from smaller radii (down the the effective resolution limit) 
would lead to even shallower inner profiles. Given our final aim to compare the simulated DM 
profiles with the GeV excess data (that require an inner slope as steep as 1.26), our choice 
is truly “conservative” in a two-fold way: (a) the maximal asymptotic slope extrapolation 
guarantees that, at any radius r, profiles steeper than the maximal asymptotic slope are not 
allowed by the simulation data; (b) the choice of the Power radius guarantees that profiles 
steeper than the maximal asymptotic slope at the Power radius are not allowed by the 
simulation data within the resolution/convergence limit of the simulation. 

The profile so determined is the steepest power-law that can be accommodated by the 
simulation within its resolution/convergence limit. 

The maximal asymptotic slope values in the simulation we analysed are: 

• EAGLE IR (10 haloes): 0.89 < 7max < 1-95, at Rpo 3 = 3.6 kpc 

• EAGLE HR (2 haloes): 0.94 < 7max < 0.98 at Rpo 3 = 1-8 kpc 

• APOSTLE IR (2 haloes): 0.50 < 7max < 0.62 at Rpo 3 = 1-8 kpc. 

Only among the relatively low-resolution EAGLE IR galaxies we find objects with a 
maximal slope as steep as the one required to explain the GeV excess emission, 1.26 ± 0.15, 
and that is only because the maximisation of the slope, if calculated from the relatively large 
Power radius of 3.6 kpc, leaves considerable freedom in the choice of the profile. For all 
other DM haloes (of EAGLE IR, EAGLE HR and APOSTLE IR selected MW-like galaxies) 
the maximal asymptotic slopes at the Power radii are significantly shallower than the steep 
profile required to fit the data. The result of our extrapolation indicates that no simulated 
halo has enough DM mass within the Power radius to support profiles as steep as 
(as required by the data). In the following discussion, we will consider only high resolution 
haloes, namely EAGLE HR and APOSTLE IR, since we believe that for those even a very 
conservative approach might lead to realistic results, contrary to what happens for EAGLE 
IR, where the resolution is too low to obtain useful constraints. 

We adopt the spherically averaged (and extrapolated) DM density profiles of the 4 
haloes in EAGLE HR and APOSTLE IR to perform a joint fit in the ten regions analysed 
by [25]. The aim is to assess how well the GeV excess data are compatible with the spherically 
symmetric DM signal predicted by EAGLE simulations. The DM signal is the one predicted 
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by the DM density distribution of the 4 selected MW-like galaxies, adopting the maximal 
asymptotic slope extrapolation. For a generic WIMP DM candidate, the predicted flux of 
photons is written as: 


He 


(av) dN^ 


Stt dE 


/ ds p‘^{r{s,'ilj)). 

' l.o.s. 


(5.2) 


The particle physics parameters {{crv) and mass) are determined by the fitting procedure (see 
below). The integral along the line of sight (l.o.s.) of the DM density profile squared (which 
depends on the l.o.s. variable s and on the angle from the GC, ^p), the so-called J-value, in¬ 
stead is fixed for each selected galaxy. We fix the branching ratio to be 100% annihilation into 
b-quarks (our conclusions are however independent of the assumed annihilation channel) and 
we compute the corresponding gamma-ray spectrum, dN^/dE, with DarkSUSY 5.1.1 [80]. 
The fitting procedure fully takes into account the background model systematic uncertain¬ 
ties in each of the ten sub-regions, but it does not account for possible correlations between 
different sub-regions. The inclusion of those correlated uncertainties has been demonstrated 
to allow more freedom for models fitting the excess, as thoroughly shown in the case of DM 
in ref. [27]. The function is; 


10 24 

= YjYj - Pik) , (5.3) 

2=1 j^k=l 


where dij (pij) represents the measured (predicted) flux in the segmented region i and ener¬ 
getic bin j ; S* ^ is the covariance matrix for energy bins j and k in region i [25] . The p-values 
quoted below refer to this function for the ten sub-regions fit, assumed to follow a xt 
distribution with k = 240 — 2 degrees of freedom. 

The free parameters of the fit, {av) and DM mass, are determined by minimising the x^ 
in eq. 5.3. We emphasise that any difference in the absolute normalisation of the DM density 
profile is re-absorbed in the fitting procedure by a different combination of annihilation cross- 
section and DM mass that still provides a good fit to the data. The goodness of fit is only 
determined by the shape of the profile. 

In table 3, we quote the best-fit values for (av) and mass, as well as the p-values, for 
the best-fit galaxies in the EAGLE HR and APOSTLE IR run (and for a gNFW profile with 
7=1.26). Correspondingly, we show constraints in the {{av), m^) plane from the fit to the 
ten sub-regions in figure 8. As expected, owing to the lower J-value due to the flattening 
in the inner kiloparsecs, the best-fit cross section in the case of a DM density profile drawn 
from the EAGLE HR and APOSTLE IR selected MW-like galaxies is slightly higher than 
the best-fit cross section obtained when adopting a gNFW (7 = 1.26) in the fit.^ Although 
it has been shown that current upper limits on the annihilation cross section from dwarf 
spheroidal galaxies are not yet in tension with the DM parameter space preferred by the DM 
interpretation of the GeV excess when a gNFW profile with 7 = 1.26 is assumed [27], we 
note that the results obtained by using a shallower DM profile might lead sooner to a tension 
with dwarf limits. 

Figures 9 and 10 show the fluxes in the ten sub-regions for the density profiles of the 4 
MW-like galaxies and particle physics parameters as obtained by the fit for each galaxy. In 
general, all DM haloes of selected galaxies fail in accounting for the flux in the two innermost 

®We note that an additional source of discrepancy in the overall normalisation of the contour regions in 
figure 8 is due to the different local dark matter density values. 
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regions, that extend up to 5° (0.75 kpc) above and below the Galactic plane. This is due to 
the flattening of the profile on scales of 1 - 2 kiloparsecs, even when extrapolating the density 
profile with the maximal asymptotic slope method. On the other hand, in the outermost 
regions the fit is relatively good because the profile preferred by the data and that used by 
the model are similar. 

In figure 11 we display the angular profile of the DM signal as predicted when using the 
density profile of the best-fit MW-like galaxies in the EAGLE HR and APOSTLE IR runs. 
The set of 4 lines refers to the 4 galaxies. Eor galaxies in the same simulation run the particle 
physics factor is fixed to the best-fit parameters of the best-fit galaxy of that simulation. The 
scatter between two lines of the same colour thus indicates the uncertainty affecting the J- 
value, including the difference in the local DM density. Eor comparison, the expected DM 
signal for the gNEW profile with 7=1.26 is also shown. The data point represents the flux 
of the GeV excess at 5° from the GC. It was demonstrated in [25] that this is a good pivot 
point since the flux only mildly depends on the inner slope of the DM profile. Below 5° 
the expected signal from EAGLE HR and APOSTLE IR selected MW-like galaxies does not 
account for the extra-emission at the GG. The mild discrepancy between all simulated haloes’ 
angular profiles and the gNEW prediction for R > 2.4 kpc is due to the different local density 
normalisation. 

Given the DM density profile predicted by the hydrodynamic simulation (and the fact 
that those profiles are not well described by any of the parameterisations commonly adopted 
in the literature), we re-analyse Fermi-LAT gamma-ray data using a GeV excess template 
built on the spatial distribution of the galaxy that best fits the GeV excess data, i.e. the 
best-fit galaxy of the EAGLE HR run. In particular, we are interested in a quantitative 
estimate of the goodness of fit of such a DM template with respect to the best-fit template 
adopted in [25]. As thoroughly explained in [25], such an analysis depends on the Galactic 
diffuse model adopted. We run the fitting template procedure described in [25] for the DM 
template based on the EAGLE HR DM halo profile extrapolated with maximal asymptotic 
slope at the Power radius of 0.94. We analyse all 60 models for the Galactic diffuse emission 
presented in [25]. We find that the Galactic diffuse emission model leading to the best fit 
for the EAGLE HR DM halo template is the so-called model E, the same foreground model 
giving the best-fit result for a gNEW template with 7 = 1.26. In figure 12, we show the 
variation of the test statistic® TS(= —21n£), ATS, as function of the energy for the best-fit 
Galactic diffuse model for the EAGLE HR DM halo template when compared to the best-fit 
one for the gNEW template with 7 = 1.26. The ATS is an indicator of how much better 
a model performs with respect to another at all energies above 1 GeV. It is evident that in 
general the EAGLE HR DM template leads to a worse fit in the energy range relevant for 
the GeV excess (1-3 GeV), while it performs similarly to the gNEW template at higher 
energy. To give a term of comparison for the ATS values, we show in the same figure 
the ATS between the next-to-best-fit and the best-fit Galactic diffuse emission model when 
using the gNEW template with 7 = 1.26. The variation in ATS due to the different Galactic 
diffuse emission model are as large as 50, although the differences are more pronounced at 
low energies. Remarkably, the uncertainties due to the variation of the DM template are 
comparable to the ones due to the Galactic diffuse modelling. 

®The TS is a measure of the goodness of fit, being directly related to the likelihood. 
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6 Conclusions 


In this work, we have used the set of cosmological, hydrodynamic simulations of the EA¬ 
GLE [41] and APOSTLE [49] projects. We extracted the simulated haloes in the generous 
Milky Way-like mass range of 5 x 10^^ < M200/MQ < 1 x 10^^. Eor each resolution run, we 
required that the galaxies satisfy observational properties of the Milky Way, besides the un¬ 
certain constraint on the halo mass. In particular, we defined galaxies to be Milky Way-like 
if (i) they give a good ht to the Galactic rotation curve [5], (ii) have a stellar mass in the 
3cr observed Milky Way stellar mass range [53], and (hi) show a dominant disc in the stellar 
component. We then analysed the dark matter density profiles of the selected Milky Way-like 
galaxies and we derived implications for the dark matter interpretation of the Fermi GeV 
excess. We summarise our findings below: 

(i) The adopted selection criteria proved to be very powerful in reducing the large variation 
in the rotation curves predicted by simulated galaxies selected on the basis of the virial 
mass only (cf. section 3). 

(ii) As a consequence, the dark matter density profiles of the Milky Way-like galaxies in 
our final selection are remarkably consistent with each other, and they have local dark 
matter densities easily compatible with the measured value of about 0.4 GeV/cm^. 

(hi) The subgrid physics model implemented in EAGLE (and APOSTLE) simulations leads 
to the dark matter density prohles of Milky Way-like galaxies exhibit a flattening in 
the inner 1.5 - 2 kpc (probably due to starburst events at low redshift) and a regime 
of contraction, from 1.5-2 kpc up to about 10 kpc, where the profile’s slope is steeper 
than 1. While the flattening in the EAGLE IR simulation might still be affected by 
numerical effects, the flattening seen in EAGLE HR and APOSTLE IR appears to have 
a physical origin, cf. section 4. 

(iv) We use the dark matter density profiles of the selected Milky Way-like galaxies to 
compute the predicted gamma-ray flux from dark matter annihilation. We extrapolate 
the dark matter prohles down to 0.3 kpc through the maximal asymptotic slope method. 
Even under the very conservative assumption of extrapolating from the Power radius, 
when performing the ht to the GeV excess data [25], we found that those dark matter 
prohles fail to reproduce the right morphology of the excess in the innermost regions 
(within 5° above and below the Galactic plane).' On the other hand, all selected Milky 
Way-like galaxies give a good ht to the GeV excess data in the other regions of interest. 
Moreover, the parameter space for annihilation cross section and dark matter mass are 
well in agreement with previous hndings within 20%. 

We thus showed that the dark matter density prohles of our selected Milky Way-like 
galaxies lead to a gamma-ray hux from dark matter annihilation in the main Galactic halo 
that cannot fully account for the observed Fermi GeV excess. Our results do not exclude 
the possibility that dark matter annihilations contribute to the GeV excess, but if that is the 
case, they require an additional component of emission in the innermost few degrees, e.g. in 
the form of point-like sources whose huxes are too dim to be detected by the Fermi LAT 
telescope as individual objects, as recently proposed in Refs. [38, 39]. 

^However, we warn the reader that, at the moment, we cannot really rule out that the shallower profile 
towards the centre is due to numerical or subgrid physics effects [41]. 
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Figure 9. Fluxes in the ten sub-regions for the two MW-like galaxies selected for the EAGLE HR 
run. The solid line indicates the best-fit galaxy, while the dashed line corresponds to the other galaxy 
in the final sample (in most of the regions the two lines overlap because of the fit results). For each 
galaxy, the flux is computed for mass and cross section as obtained by the fit to the ten sub-regions 
(see text for details). The black points are the GeV excess data as derived in [25], with corresponding 
systematic uncertainties {grey boxes). The insets refer to the region of interest over which the flux is 
averaged, see [25] for the definition of the ten regions. 
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Figure 10. Same as figure 9 for the two MW-like selected haloes in the APOSTLE IR run. 
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Figure 11. The magenta (orange) lines represent the angular profile of the DM signal predicted by 
the selected MW-like halo’s density profiles in the EAGLE HR (APOSTLE IR) run. The solid line 
refers to the best-ht halo, while the dashed line corresponds to the density profile of the other halo 
selected normalised to the best-fit cross section and mass of the best-fit halo. The difference between 
two lines of the same colour thus indicates the uncertainty affecting the J-value. The black dotted line 
represents the expectation for the gNFW profile with 7=1.26 - that best hts the GeV excess data 
while the black point is the flux of the GeV excess at 5° away from the GG (that was shown to 
be a good pivot point in [25]). The vertical grey dashed line represents the Power radius for the two 
simulation runs. 
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Figure 12. Variation of the TS(= —21n£) between the best-fit Galactic diffuse emission model for 
the EAGLE HR DM profile template and the best-fit Galactic diffuse emission model for the gNFW 
template [25]: the blue solid line uses as DM template the EAGLE HR density profile extrapolated 
with the maximal asymptotic slope method from the Power radius down to smaller scales. The TS 
variation of the next-to-best-fit Galactic diffuse emission model for the gNFW template with respect 
to the best-fit model is shown by the red line. 
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Figure 13. DM density profiles for the two selected haloes with the correct total stellar mass (within 
3 (t) and with the lowest and largest halo mass, M 200 , for the EAGLE IR {left panel) and EAGLE HR 
{right panel) runs. 

A Variation with halo mass 

In figure 13, we display the DM density profile for two haloes in the EAGLE IR and EAGLE 
HR runs that have the correct MW observed stellar mass range (within 3cr) and that have 
the lowest and highest M 200 respectively. This figure illustrates that the features of the DM 
density profile (i.e. the inner flattening and the contraction regime) do not strongly depend 
on M 200 - While for EAGLE IR also the normalisation of the density profiles is well matched, 
for EAGLE HR haloes the overall profile normalisation might differ up to a factor of 5 at 
the resolution limit. Nevertheless, for the purposes of our analysis, such a normalisation 
discrepancy is not critical. Indeed, as explained in section. 5, the goodness of fit to the 
GeV excess data is driven by the shape of the DM density profile in the range 0.3 - 3 kpc 
and by the extrapolated maximal asymptotic slope. To demonstrate that the halo mass 
variation does not introduce any bias in our analysis, we perform the fit to the GeV excess 
data for the halo with the largest mass in figure 13 (right panel). The best-fit parameters are 
(av) = (2.45 ± 0.20) x 10“^®cm^/s and = 45.4 ± 2.1 GeV. Those values are compatible 
(within 2a) with the EAGLE HR entry in table 3. In light of this result, we are confident 
that discrepancies in the halo mass does not alter our final conclusion. 

B Varying local Galactic parameters 

Throughout this paper, we have assumed Rq = 8 kpc, vq = 230 km/s, and the peculiar 
velocity of the Sun, {U,V,W)q = (11.10,12.24,7.25) km/s given in Galactic coordinates. 
In this appendix, we consider different choices of Rq, vq and Vq, showing that the main 
conclusions of the paper remain unchanged. We follow ref. [81], and consider the following 
four different configurations: 

(a) Rq = 7.98 kpc, vq = 214.44 km/s for Vq = 26 km/s; 

(b) Rq = 7.98 kpc, vq = 236.94 km/s for Vq = 5.25 km/s; 

(c) Rq = 8.68 kpc, Vq = 235.53 km/s for Vq = 26 km/s; 

(d) Rq = 8.68 kpc, vq = 258.19 km/s for Vq = 5.25 km/s. 
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(a) 



(c) 



(b) 



(d) 



Figure 14. Observed MW angular circular velocity as a function of galactocentric distance from [5] 
{black points and error bars) for configuration (a) {top left), (b) {top right), (c) {bottom left), and (d) 
{bottom right), as well as the best-fit halo in the EAGLE HR {magenta) and EAGLE IR {green) runs. 


Configuration 

xV(A^- 1), EAGLE HR 

xV(A^- 1), EAGLE IR 

(a) 

70.48 

15.88 

(b) 

73.21 

16.15 

(c) 

98.49 

16.83 

(d) 

93.62 

17.92 


Table 4. Reduced values for the best ht haloes which satisfy the criteria in section 3 in the 
EAGLE HR and EAGLE IR runs for the four configurations discussed in the text. 

These configurations are used to obtain the observed MW rotation curve data. We 
perform the analysis discussed in section 3.1 to find the best-fit haloes which satisfy the 
criteria in section 3. The observational data for the MW angular circular velocity as well as 
the best-fit haloes for the EAGLE HR and EAGLE IR runs are shown in figure 14. The four 
panels of the figure correspond to the four configurations. The reduced x^ for these haloes 
are specified in table 4. The DM density profiles of the MW-like selected haloes are shown 
in figures 15 and 16 for the EAGLE IR and EAGLE HR runs, respectively. All haloes (both 
in the EAGLE IR and EAGLE HR runs) show a significant deviation from the NEW profile 
and the common features of the density prohle discussed in section 4 can be found also when 
varying Rq, vq and V© parameters. 
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Figure 15. Same as figure 7 for the EAGLE IR run and for the four combinations of vg and Bg. 
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Figure 16. Same as figure 7 for the EAGLE HR 



run and for the four combinations of vq and Rq. 
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